A Guide for Diagnostics in Linear Mixed-Effects Models
Introduction to This Page
Linear mixed-effects models (LMM) are a generalization of the classical linear regression models to account for grouping effects. They are particularly useful in settings where repeated measurements are made on the same statistical units, such as in a longitudinal study. LMMs are popular in a wide variety of practical disciplines for various kinds of inference purposes. Due to the flexibility and complexity of LMMs, the diagnostics techniques are sloppily documented. Therefore, this page collects in one place the procedures and various tools for model diagnostics in linear mixed-effects models, along with discussions and analysis from a statistical perspective.
Model Specificiation
The standard representation in univariate form for a linear mixed-effects model (LMM) is \[ \begin{aligned} &Y_j=X_j\beta+Z_j b_j+\varepsilon_j\\ &b_j\perp \varepsilon_j,\quad b_j\sim N_q(0,D),\quad \varepsilon_j\sim N_{n_j}(0,R_j), \end{aligned} \] where \(j=1,...,J\) is the index of the groups, and \(n_j\) is the sample size for group \(j\).
Model Assumption and Diagnostics
Including random effects makes a standard linear regression model becomes a linear mixed-effects model. Then, similarly to standard regression models, the model diagnostics for linear mixed-effects models mainly focus on the 4 assumptions: 1) linearity, 2) normality, 3) constant variance, and 4) influence measure and outliers. As will be seen in the following section, the diagnostics tools are similar to those in the standard regression model.
What makes the diagnostics for linear mixed-effects models special?
Flexible model specifications and unstructured covariance matrix
The performance of the classical linear regression models highly relies on the conformity to the seven classical assumptions which make the OLS estimates good. However, for linear mixed-effects models, the model specification can be highly flexible.
More specifically, for the error term \(\varepsilon_j\sim N_{n_j}(0,R_j)\), although we typically simplify the model with the constant variance assumption, i.e., \[ \mathbb{V}ar(\varepsilon_j)=\sigma^2 I_{n_j}, \] an elaborate structure for \(R_j\) can also be desirable when, for example, \(R_j\) depends on the value of a covariate, or when it is longitudinal study so that an autoregressive structure is needed [1]. Then, the diagnostics should be adjusted accordingly.
Residuals within each group
For linear mixed-effects models, two steps can be considered in diagnostics: one is for the residuals for all the observations as a whole, and the other is for the residuals for each of the groups, to see if any of the groups show different patterns and violate the model assumptions. Particularly, as mentioned above, if the distributional specifications are different for different groups, then the diagnostics for residuals as a whole are not valid; but instead, we should analyze the residuals patterns for each group.
However, sample size issues can be further raised when analyzing the group-specific residuals. That is, we may have some groups with small sample sizes (and in fact, this can often be expected to see in LMMs because one of the key reasons to include random effects is to apply shrinkage and to borrow information to build more reliable estimates for those groups with few samples), and the diagnostics for the true unobservable errors are difficult in small samples by examination of residuals [2]. Therefore, to have applicable and reliable diagnostics tools which enable us to check the appropriateness of our model, and also to avoid running into the identifiability issues, typically we tend to consider a simpler model specification as long as it is still reasonable, except when we have strong background knowledge or belief.
Higher-level variables and lower-level variables
Linear mixed-effects models estimate the effects of both lower-level and higher-level variables. The influence diagnostics is different for linear mixed-effects models in that there are multi-level influences involved. It is possible that in some cases a higher-level group is influential, while in others, a single observation within a group is influential. Thus, the investigation into both the higher- and the lower- influences is necessitated [1]. In addition, in linear mixed-effects models, the influence measure can be of particular interest to identify groups that are different from others. See the section below for more detailed discussions on the influence diagnostics in linear mixed-effects models.
Simulation Setup
For illustration purpose, in this report we consider one of the simplest examples of linear mixed-effects models given by \[ \begin{aligned} &Y_{ij}=(\beta_0+b_{0,j})+(\beta_1+b_{1,j})X_{ij}+\varepsilon_{ij},\\ &b_{0,j}\sim N(0,\tau_0^2)\perp b_{1,j}\sim N(0,\tau_1^2)\perp \varepsilon_{ij}\sim N(0,\sigma^2),\quad i=1,...,n_j,\; \;j=1,...,J. \end{aligned} \] That is, the linear mixed-effects model considered in this report includes random intercepts, and random slopes for a single predictor \(X\), as well as a fixed overall intercept and a fixed overall slope.
In the most basic setting, we will generate \(J=9\) groups, where 3 contain 5 observations, 3 contain 15 observations, and the last 3 contain 150 observations. Besides, the parameters are chosen to be \(\beta_0=0,\beta_1=1,\tau_0^2=5,\tau_1^2=10,\sigma^2=5\).
Linearity
Linearity assumption requires that the relationship between the response and the predictors is linear. In other words, we have a mean function that is linear in the coefficients, i.e., \[ \mathbb{E}(Y_i|X_i,Z_i)=X_i\beta+Z_i b_i. \] Theoretically, if the mean function (and other assumptions) are correct, then all possible plots of residuals versus any linear function of the predictors should resemble a null plot [2]. That is, the linear predictor has captured the information in the response, leaving nothing but the random noise (which is assumed to be normally distributed).
Since the residuals are orthogonal to the column space of the design matrix, then all linear combinations of the predictors can be plotted against the residuals to check the linearity. Nonetheless, the most commonly used diagnostic tool is the plot of residuals versus fitted values, while sometimes we may also refer to the plots of residuals versus each of the predictors. Besides, as mentioned above, for linear mixed-effects models, we can look into the residual plots both with residuals mixed together and with residuals grouped by the grouping variables to see if the linearity assumption is not met in any particular groups.
As shown in the above plots, the plot of the residuals versus the fitted values and the plot of the residuals versus the predictor \(X\) both resemble null plots, confirming that the linearity assumption is met. This is a natural result under our simple and ideal simulation data. Also, in our example we only include a single predictor \(X\), so that the plot against the \(X\) and the plot against the fitted values are highly similar.
Notably, in the above plots, the groups with small sample sizes (i.e., “g1”, “g2”, and “g3”) and even the groups with moderate sample sizes (i.e., “g4”, “g5”, and “g6”) don’t really resemble null plots. As mentioned above, this is the case where the residual plots with only few scatters are not quite informative and reliable. Under such circumstances, we may refer to the residuals for the larger groups such as “g7”, “g8”, and “g9”, or refer to the residuals mixed together to draw our conclusions.
Normality
The normality assumption requires that the error term follows a normal distribution, i.e., \[ \varepsilon_j\sim N_{n_j}(0,R_j) \] The most commonly used diagnostics tool for normality is the normal probability plots (aka the QQ plot).
As shown in the above plot, again, we can look into the QQ plot for all the residuals and the QQ plot for residuals in each of the groups. Scatters following a 45° straight line indicates the conformity to the normality assumption.
Constant Variance
The constant variance assumption (aka homoscedasticity) requires that the error terms have the same variance across the data. Note that linear mixed-effects models do not necessarily assume constant variance but often favor a more elaborate structure for the error variance, as mentioned above; for example, autoregression structure is desirable for longitudinal study. The following diagnostics is valid under the simplified model assumption \(\varepsilon_{ij}\stackrel{iid}{\sim} N(0,\sigma^2)\), where plots of standard residuals versus fitted values or versus the predictors can be used for diagnostics.
As shown above, the pattern that the variation of the standardized residuals remains constant across the fitted values or the predictor values indicates that the constant variance assumption has been met. Again, similar to what has been seen in the linearity diagnostics, groups with small sample sizes are not appropriate to draw conclusion on the model assumption.
In addition, for linear mixed-effects models, the box plot of standardized residuals by groups is another useful tool helping to quickly identify groups with specially high or low variation and compare the distribution of the standardized residuals across groups, as shown in the above plot.
Influence
Influence measure
Single cases or small groups of cases can strongly influence the fit of the model. Cases whose removal causes major changes in the analysis are called influential. Basically, we require that all data points have the same amount of influence on our parameter estimates, with a hope that the model conclusions can be robust to slight perturbations or missingness of part of the data. Based on idea of deletion diagnostics [1], we’ve built several influence diagnostics tools, among which the two commonly used in LMMs are the Cook’s distance and DFBETAS.
DFBETAS:
DFBETAS (standardized difference of the beta) measures the level of influence observations have on single parameter estimates. It is calculated as the difference in magnitude of the parameter estimate between the model including and the model excluding the group, standardized by dividing by the standard error of the estimate that excludes the group (to prevent variance inflation from masking the level of influence). Formally, the DFBETAS for group \(j\) and parameter \(k\) is defined as \[ DFBETAS_{jk}=\frac{\hat{\gamma}_k-\hat{\gamma}_{k(-j)}}{se(\hat{\gamma}_{k(-j)})}, \] where \(\hat{\gamma}_k\) is the original estimate of the \(k\)-th parameter, and \(\hat{\gamma}_{k(-j)}\) is the estimate of the same parameter after group \(j\) has been excluded from the data [1].
A commonly used cutoff for identifying overly influential observations is \(\frac{2}{\sqrt{n}}\), where \(n\) is the number of groups at the level of removal \((-j)\) [1]. As shown in the plot below, (the absolute value of) the estimates for each coefficients are compared with the cutoff. It turns out that group “g7” strongly influences the estimation of the slope of predictor \(X\).
Influence diagnostics with DFBETAS will become a difficult task when the number of parameters is large. When we compute the DFBETAS for each covariate \(k\) for the group \(j\), it is possible that the group have strong influence on some covariates while not influential on the others. Then it is hard to decide on the overall influence of that group. Besides, if the number of groups and the number of parameters are both large, then we will have too many DFBETAS to wade through. Typically, DFBETAS is useful when the number of parameters is not quite large.
Cook’s distance:
Cook’s distance gives us a summary for influence on all parameter estimates, which is a better choice when the number of parameters is large. Formally, it is defined as \[ C_j=\frac{1}{p}(\hat{\gamma}-\hat{\gamma}_{(-j)})^T\hat{\Sigma}_{(-j)}^{-1}(\hat{\gamma}-\hat{\gamma}_{(-j)}), \] where \(p\) is the length of \(\beta\), and \(\hat{\Sigma}_{(-j)}\) is the covariance matrix of the parameter estimates excluding group \(j\). The commonly used cutoff is \(\frac{4}{n}\) where \(n\) is the number of groups in the grouping factor being evaluated. If there is just one parameter in the model, then Cook’s distance is the square of DFBETAS for that parameter [1].
The above plots show the Cook’s distance measure at the group level. Groups with Cook’s distance greater than the threshold are identified as influential. Then, the plot of the fitted values versus the observed response values, or the plot of the response values versus any of the predictors, can be used to further examine how exactly the influential groups are different from others. Additionally, in the above plots (and also in the following plots), in each of the plots of the fitted values versus the observed values, a straight line with intercept 0 and slope 1 is added for comparison, while in that of the response versus the predictor, the regression line with the corresponding slope and intercept estimated by the LMM is added. Also, the sample sizes of each group are included in parenthesis as supplementary information for analysis.
Similarly, at the lower individual level, we can compare the Cook’s distance of each of the individual point to the corresponding threshold. Furthermore, we can look into the plot of the fitted values versus the observed response values, or the plot of the response values versus any of the predictors, to identify how the points identified to be influential actually influence the model estimates, as shown in the above plot.
Influence interpretation
As mentioned above, the influence diagnostics is complex for linear mixed-effects models in that both higher-level and lower-level influences should be examined. On the other hand, it is also of particular interest to conduct the influence diagnostics as the higher-level analysis helps to identify groups that are different from others, especially for the LMM that is built with a focus on the heterogeneity among groups.
Then, unlike in the standard regression model where the influence analysis mainly focuses on the individual cases, in LMMs, the interpretation for a group being influential can be complicated. In addition to analyzing how the influential groups truly differ from others in nature and whether the influential groups are nuisance outliers, the other factors worth consideration include: 1) the group size, and 2) the number of groups in the data set. For example, one possibility is that when there are not many groups, it is likely to identify influential groups. As shown in the plots below, we include 4 groups, all of which have the same either small or large sample sizes. We identify 2 out of the 4 groups as strongly influential in the sense of Cook’s distance, although indeed all data come from the same data generating model.
- Group size = 10
- Group size = 100
Another possible scenario we may expect to see in practice is that we select out groups with extremely large sample sizes as the influential groups. However, it can be hard to tell the reason without a specific investigation into the data set. Sometimes it is because the group as a whole is truly unique in the coefficient estimates, while sometimes it is because the huge group merges too many cases, many of which should have been modeled differently in nature but under the current model, are identified as influential cases and also leading the whole group to be influential.
Apparently, the multiple factors contribute in a mixed way to the influence measure in linear mixed-effects models, meaning that it is often difficult for us to conclude on the decisive factor dominating the influence magnitude. In short, no definite conclusion can be made consistently for all the influence analysis, and this is where both statistical investigation and domain knowledge should come into play.
Conclusion and Discussion
This report has summarized several basic diagnostics tools for linear mixed-effects models. To be sure, the real data going into a linear mixed-effects model can be much more complex and involved, hard to be accommodated through our simulation and illustration. Nonetheless, this document is still helpful as a detailed guide for model diagnostics for LMMs. In this final section, there are some additional points that we would like to mention.
Similar to the standard regression model, the linear mixed-effects models can be generalized to fit the non-normal response. Then, in the generalized linear mixed-effects models, the model diagnostics procedures and tools will be different accordingly.
Last but not least, as we have emphasized several times, the model specification can be highly flexible for linear mixed-effects models. Thus, in practice, for a linear mixed-effects model we may not often expect to see an exhaustive and strict diagnostics procedure which is usually desired for a standard regression model. Instead, we may choose to focus on one or more specific diagnostics, based on the our modeling purpose and the research questions of interest.
Reference
[1] Akande, O. (2021). Slides for STA610L: Multilevel and Hierarchical Models (Spring 2021), Duke University. https://sta-610l-s21.github.io/Course-Website/
[2] Weisberg, S. (2014). Applied Linear Regression, Fourth Edition. Hoboken, New Jersey: John Wiley & Sons, Inc.
[3] Loy, A., Hofmann, H., & Cook, D. (2017). Model choice and diagnostics for linear mixed-effects models using statistics on street corners. Journal of Computational and Graphical Statistics, 26(3), 478-492. https://doi.org/10.1080/10618600.2017.1330207.
Appendix
All simulation results in this report are reproducible with the codes below.
library(tidyverse)
library(lme4)
library(influence.ME)
library(ggpubr)
set.seed(532)
# color to be used
color1 <- "#474747" # dark
color2 <- "#4e89ae" # blue
color3 <- "#e46457" # red
# true parameters
beta0 <- 0 # intercept
beta1 <- 1 # slope
tau2.0 <- 5 # variance of random intercepts
tau2.1 <- 10 # variance of random slopes
sigma2 <- 5 # variance of residual
# generate data
generate_data <- function(J, J.large, J.small,
n.small = 5, n.mid = 15, n.large = 150) {
# J: number of groups
# J.small, J.large: number of small groups and large groups
# n.small, n.mid, n.large: sample size of the small, middle, and large groups
group_size <- c(rep(n.small, J.small), rep(n.mid, J-J.large-J.small),
rep(n.large, J.large)) # length=J, sample sizes of each group
b0 <- rnorm(J, 0, sqrt(tau2.0)) # length=J, random intercept for each group
b1 <- rnorm(J, 5, sqrt(tau2.1)) # length=J, random slopes for each group
# b1 <- runif(J, 1, 10)
# expand the random effects , repeat each by n_j times
b0.expand <- map2(b0, group_size, ~rep(.x, times = .y)) %>% unlist()
b1.expand <- map2(b1, group_size, ~rep(.x, times = .y)) %>% unlist()
labelpart <- paste0("g", seq_len(J))
if (J >= 10) {
for (i in seq_len(9)) {
labelpart[i] <- paste0("g0", i)
}
}
group_label <- map2(labelpart, group_size, ~rep(.x, times = .y)) %>% unlist()
epsilon <- rnorm(sum(group_size), 0, sqrt(sigma2))
X <- runif(sum(group_size), -5, 5)
Y <- (beta0 + b0.expand) + (beta1 + b1.expand) * X + epsilon
# add perturbation to mimic the real data
rand_index <- sample(sum(group_size), ceiling(sum(group_size) / 5))
Y[rand_index] <- Y[rand_index] + rnorm(length(rand_index), 0, 3)
data.frame(Y, X, group = group_label)
}
# modeling
get_model <- function(data) {
lmer(Y~ (1+X|group) + X, data)
}
# DIAG1: influential
# group-level cood D
get_cook.group <- function(model) {
model.inf.group <- influence(model, group="group")
cook.group <- data.frame(
cook = cooks.distance(model.inf.group),
cutoff = 4/length(unique(data$group))
) %>%
arrange(desc(cook)) %>%
mutate(rank = seq_len(nrow(.)),
influential = cook > cutoff) %>%
rownames_to_column("group")
cook.group
}
# individual-level cook D
get_cook.ind <- function(model) {
model.inf.ind <- influence(model, obs = T)
cook.ind <- data.frame(
group = data$group,
cook = cooks.distance(model.inf.ind),
cutoff = 4/nrow(data)
) %>%
mutate(influential = cook > cutoff)
cook.ind
}
# combine the original data with the cook D
get_data.withCD <- function(cook.group, cook.ind, data, model) {
# get coefficients estimates
# (`coefficients()` return the estimates of fix + random)
result.coefall <- coefficients(model)$group %>%
rownames_to_column("group") %>% rename(slope = X)
data.withCD <-
# join the group-level cook D with the original data set
inner_join(data, cook.group, by = "group") %>%
# get the fitted values
mutate(Yfitted = fitted(model)) %>%
# get the group sizes
inner_join(data %>% count(group), by = "group") %>%
# make label for plots
mutate(label = str_c(group, " (", n, ")")) %>%
# combine the coefficients estimates
inner_join(result.coefall, by = "group") %>%
# get the individual-level cook D (rename, g for group, i for individual)
rename(cook.g = cook, cutoff.g = cutoff, influential.g = influential) %>%
cbind(cook.ind[-1]) %>%
rename(cook.i = cook, cutoff.i = cutoff, influential.i = influential)
data.withCD
}
# plot: influence diagnostics
plot.inf.group <- function(cook.group, cook.ind) {
# group-level cook D
cook.group %>%
ggplot(mapping = aes(color = influential)) +
theme_classic() +
geom_hline(aes(yintercept = cutoff),
linetype = "dashed", size = 0.7, alpha = 0.4) +
geom_segment(aes(x = reorder(group, -rank), xend = reorder(group, -rank),
y = 0, yend = cook), alpha = 0.7) + coord_flip() +
geom_point(aes(x = reorder(group, -rank), y = cook), alpha = 1, size = 1) +
scale_color_manual(values = c(color1, color2)) +
labs(y = "Group-level cook distance", x = "Group",
color = "Influential Group") +
theme(legend.position = "top")
}
plot.inf.group.scatter <- function(data.withCD, y_axis, x_axis) {
# scatter plots to analyze the group-level cook D
# choice of (y_axis, x_axis): ("Yfitted", "Y") or ("Y", "X")
if (y_axis == "Yfitted" & x_axis == "Y") {
p <- data.withCD %>%
ggplot() +
geom_abline(slope = 1, intercept = 0, alpha = 0.5)
}
if (y_axis == "Y" & x_axis == "X") {
p <- data.withCD %>%
ggplot() +
geom_abline(aes(slope = slope, intercept = `(Intercept)`), alpha = 0.5)
}
p +
geom_point(aes(x = get(x_axis), y = get(y_axis), color = influential.g),
size = 0.8, alpha = 0.7) +
scale_color_manual(values = c(color1, color2)) +
facet_wrap(~label) + theme_classic() +
labs(color = "Influential Group", x = x_axis, y = y_axis)
}
plot.inf.group2 <- function(cook.group, cook.ind, data.withCD) {
# arrange all plots for group-level Cook D
p1 <- plot.inf.group(cook.group, cook.ind) + theme_classic(base_size = 12)
p2 <- plot.inf.group.scatter(data.withCD, "Yfitted", "Y") +
theme_classic(base_size = 10) + theme(legend.position = "none")
p3 <- plot.inf.group.scatter(data.withCD, "Y", "X") +
theme_classic(base_size = 10) + theme(legend.position = "none")
ggarrange(p1, ggarrange(p2, p3, nrow = 1),
nrow = 2, common.legend = T, legend = "top")
}
plot.inf.ind <- function(data.withCD, y_axis, x_axis) {
# individual-level cook D, within group scatter plot of y_axis vs x_axis
# choice of (y_axis, x_axis): ("Yfitted", "Y") or ("Y", "X")
groups_withInfInd <- data.withCD %>% count(group, influential.i) %>%
filter(influential.i) %>% pull(group) # groups with influential individuals
if (y_axis == "Yfitted" & x_axis == "Y") {
p <- data.withCD %>% filter(group %in% groups_withInfInd) %>%
ggplot(mapping = aes(color = influential.i, size = influential.i)) +
geom_abline(slope = 1, intercept = 0, alpha = 0.5)
}
if (y_axis == "Y" & x_axis == "X") {
p <- data.withCD %>% filter(group %in% groups_withInfInd) %>%
ggplot(mapping = aes(color = influential.i, size = influential.i)) +
geom_abline(aes(slope = slope, intercept = `(Intercept)`), alpha = 0.5)
}
p +
geom_point(aes(get(x_axis), get(y_axis), size = influential.i)) +
scale_color_manual(values = c(color1, color2)) +
scale_size_manual(values = c(1, 1.5)) +
theme_classic(base_size = 11) +
labs(size = "Influential Point", color = "Influential Point",
x = x_axis, y = y_axis) +
facet_wrap(~group) + theme(legend.position = "top")
}
plot.inf.ind2 <- function(data.withCD) {
p4 <- plot.inf.ind(data.withCD, "Yfitted", "Y")
p5 <- plot.inf.ind(data.withCD, "Y", "X")
ggarrange(p4, p5, nrow = 2, common.legend = T, legend = "top")
}
plot.inf.debetas <- function(model) {
model.inf.group <- influence(model, group="group")
dfbetas(model.inf.group) %>% as.data.frame() %>%
mutate(cutoff = 2/sqrt(nrow(.))) %>% rownames_to_column("group") %>%
pivot_longer(`(Intercept)`:X, names_to = "predictor") %>%
mutate(influential = abs(value) > cutoff) %>%
ggplot() +
geom_hline(aes(yintercept = cutoff), linetype = "dashed", alpha = 0.7) +
geom_hline(aes(yintercept = -cutoff), linetype = "dashed", alpha = 0.7) +
geom_hline(yintercept = 0, alpha = 0.3) +
geom_point(aes(x = group, y = value, color = influential)) +
scale_color_manual(values = c(color1, color2)) +
facet_wrap(~predictor, scale = "free") + theme_bw(base_size = 11) +
coord_flip() +
labs(y = "DFBETAS", x = "Group", color = "Influential Group")
}
## DIAG2: constant variance
get_data.withResid <- function(data, model) {
data.withResid <- data %>%
mutate(resid_scale = resid(model, scaled = T), # standardized residual
resid = resid(model, scaled = F), # residual
Yfitted = fitted(model)) # fitted values
data.withResid
}
# plot: constant variance
plot.cv.group <- function(data.withResid) {
# group-level standardized residuals boxplot
data.withResid %>%
ggplot() +
geom_boxplot(aes(x = group, y = resid_scale),
fill = "#3282b8", alpha = 0.6) +
theme_classic(base_size = 12) + labs(y = "Standardized Residuals") +
coord_flip()
}
plot.cv.ind <- function(data.withResid, x_axis) {
# individual-level standardized resid vs x_axis
# x_axis: "Yfitted" or "X"
data.withResid %>%
ggplot() +
geom_hline(yintercept = 0, size = 0.6, linetype = "dashed", alpha = 0.3) +
geom_point(aes(x = get(x_axis), y = resid_scale), alpha = 0.5, size = 1) +
geom_smooth(aes(x = get(x_axis), y = resid_scale), se = F,
color = color3, alpha = 0.7, size = 0.8) +
theme_classic() + labs(y = "Standardized Residuals")
}
plot.cv.ind.g <- function(data.withResid, x_axis) {
# individual-level standardized resid vs x_axis, grouped by group
# x_axis: "Yfitted" or "X"
data.withResid %>%
ggplot() +
geom_hline(yintercept = 0, size = 0.6, linetype = "dashed", alpha = 0.3) +
geom_point(aes(x = get(x_axis), y = resid_scale), alpha = 0.5, size = 1) +
geom_smooth(aes(x = get(x_axis), y = resid_scale), se = F,
color = color3, alpha = 0.7, size = 0.8) +
theme_classic() + labs(y = "Standardized Residuals") +
facet_wrap(~group, scale = "free")
}
plot.cv.ind2 <- function(data.withResid, x_axis) {
p5 <- plot.cv.ind(data.withResid, x_axis) + theme_classic(base_size = 9) +
theme(axis.title = element_blank())
p6 <- plot.cv.ind.g(data.withResid, x_axis) + theme_classic(base_size = 9) +
theme(axis.title = element_blank())
ggarrange(p5, p6) %>%
annotate_figure(
left = text_grob("Standardized Residuals", size = 13, rot = 90),
bottom = text_grob(x_axis, size = 13, rot = 0)
)
}
## DIAG3: normality
plot.norm.ind <- function(data.withResid) {
# individual-level qq plot
data.withResid %>%
ggplot(mapping = aes(sample = resid_scale)) +
geom_abline(intercept = 0, slope = 1, alpha = 0.6) +
geom_qq(distribution = qnorm, color = color1, alpha = 0.8, size = 0.8) +
theme_classic() +
labs(x = "Theoretical Quantile", y = "Standardized Residuals")
}
plot.norm.ind.g <- function(data.withResid) {
# individual-level qq plot, grouped by group
data.withResid %>%
ggplot(mapping = aes(sample = resid_scale)) +
geom_abline(intercept = 0, slope = 1, alpha = 0.6) +
geom_qq(distribution = qnorm, color = color1, alpha = 0.8, size = 0.8) +
theme_classic() +
labs(x = "Theoretical Quantile", y = "Standardized Residuals") +
facet_wrap(~group)
}
plot.norm <- function(data.withResid) {
p7 <- plot.norm.ind(data.withResid) + theme_classic(base_size = 9) +
theme(axis.title = element_blank())
p8 <- plot.norm.ind.g(data.withResid) + theme_classic(base_size = 9) +
theme(axis.title = element_blank())
ggarrange(p7, p8) %>%
annotate_figure(
left = text_grob("Standardized Residuals", size = 13, rot = 90),
bottom = text_grob("Theoretical Quantile", size = 13, rot = 0)
)
}
## DIAG4: linearity
plot.li.ind <- function(data.withResid, x_axis) {
# individual-level resid vs x_axis (Yfitted or Xpredictor)
# x_axis should be "Yfitted" or "X"
data.withResid %>%
ggplot() +
geom_hline(yintercept = 0, size = 0.6, linetype = "dashed", alpha = 0.3) +
geom_point(aes(x = get(x_axis), y = resid), alpha = 0.5, size = 1) +
geom_smooth(aes(x = get(x_axis), y = resid), se = F,
color = color3, alpha = 0.7, size = 0.8) +
theme_classic() + labs(y = "Residuals")
}
plot.li.ind.g <- function(data.withResid, x_axis) {
# individual-level resid vs x_axis (Yfitted or Xpredictor), grouped by group
# x_axis should be "Yfitted" or "X"
data.withResid %>%
ggplot() +
geom_hline(yintercept = 0, size = 0.6, linetype = "dashed", alpha = 0.3) +
geom_point(aes(x = get(x_axis), y = resid), alpha = 0.5, size = 1) +
geom_smooth(aes(x = get(x_axis), y = resid), se = F,
color = color3, alpha = 0.7, size = 0.8) +
theme_classic() + labs(y = "Residuals") +
facet_wrap(~group, scales = "free")
}
plot.li.ind2 <- function(data.withResid, x_axis) {
p9 <- plot.li.ind(data.withResid, x_axis) + theme_classic(base_size = 9) +
theme(axis.title = element_blank())
p10 <- plot.li.ind.g(data.withResid, x_axis) + theme_classic(base_size = 9) +
theme(axis.title = element_blank())
ggarrange(p9, p10) %>%
annotate_figure(
left = text_grob("Residuals", size = 13, rot = 90),
bottom = text_grob(x_axis, size = 13, rot = 0)
)
}